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Introducti on 


In this Report, we summarize research work done at 
Cambridge Hydrodynamics, Inc. under Contract No. NAS1-14907 
with NASA Langley Research Center. More detailed expositions 
of the work described here is given in the publications cited 
in the References. 

In Sec. 2, we summarize our work done to optimize the 
numerics of the SALLY computer code written by us and used 
at NASA LaRC and elsewhere to analyze the stability of laminar 
flow control wings. We have succeeded in speeding up earlier 
versions of SALLY by about a factor 3 while reducing the memory 
requirements by about a factor 2. 

In Sec. 3, we summarize work done by Dr. David J. Benney 
of CHI on the relation between temporal and spatial stability 
theory in laminar flow control problems. 

In Sec. 4, we describe new methods that we have developed 
for the solution of three-dimensional compressible flow 
stability problems by spectral methods. 

In Sec. 5, we describe new spectral methods for the solution 
of boundary layer equations. These new methods are both highly 
accurate and efficient. 

Finally, in Sec. 6, we describe our work extending the 
SALLY code to study the nonlinear, non-parallel stability 
theory of boundary layer flows. 
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2 . Optimal Numerics for the SALLY Stability Analysis Code 

The SALLY Computer code was developed for NASA Langley 
Research Center by CHI under contracts L-47262-A and NAS1-14427. 
It is essentially a 'black-box 1 stability analyzer for three- 
dimensional boundary layer flows over LFC wings. In the present 
work, we have speeded up the code by about a factor 3 while 
reducing core requirements by about a factor two. In addition, 
several new features and options have been added for the 
convenience of users. A detailed exposition of the numerical 
methods and results obtained by SALLY will be given in a 
forthcoming paper. Here we summarize the new developments. 

(i) Fast Eigenvalue Solvers — A new fast eigenvalue 
solver has been implemented successfully in SALLY. We use a 
cubically convergent inverse Rayleigh iteration^ to improve 
an approximate eigenvalue to within roundoff accuracy in 
about 3 iterations. Practice has shown that even a relatively 
crude guess for the eigenvalue gives good results with this 
local iteration scheme. 

The heart of this new eigenvalue solver is new assembly 
language matrix manipulation routines. Assembly language 
fully-pivoted LU decomposition and matrix multiply programs 
have been written in COMPASS for the CDC 174-176 series machines. 
These codes are now operational at NASA LaRC. Because of the 
difference between CPU speed and memory access time balances 
for these machines, it was found to be essential to write 
separate routines for each machine. For example, the Cyber 176 
has a fast CPU and a fast memory cycle allowing code optimiza- 
tion using stackloop programming techniques. On the other hand, 
the Cyber 175 has a relatively slow memory cycle so the codes 
for this machine were specially designed to minimize the memory- 
access penalty. Thus, the Cyber 175 code runs about twice as 
fast on the 175 as the 176 code runs on the 175. 
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(ii) New Global Eigenvalue Routine — The global eigenvalue 
routine of the original SALLY code^ has been replaced by an 
improved code that both runs faster and requires less memory. 

This global eigenvalue routine is essential in flow regions 
where a good guess is not available, like near the leading 

edge of the airfoil or in the stable region of the LFC wing 
where practice has shown that the character of the least 
stable mode changes greatly. 

(iii) New Group Velocity Code — A new, extremely efficient 
code for the calculation of the group velocity has been imple- 
mented in SALLY. In the original SALLY code, the group velocity 
was computed in two ways, at the option of the user. A first- 
order finite difference approximation to the group velocity was 
obtainable by computing the eigenvalue at two nearby wave-vectors; 
alternatively, an analytic computation of the eigenvalue using 

2 

the adjoint eigenfunction (determined numerically) was available. 

The new code is based in the interesting observation that the 
adjoint eigenfunction used in the analytic formula for the 
group velocity need not be the adjoint eigenfunction of the 
Orr-Sommerf eld stability equation, but rather can be the adjoint 
eigenvector of any matrix that is similar to the matrix representa- 
tion of the Orr-Sommerf eld equation. 

We have implemented the latter idea by use of the adjoint 
vectors employed in the local Rayleigh iteration scheme. In this 
case, the group velocity is immediately available upon simple 
algebraic manipulation of the eigenvector and the matrix adjoint. 
The new procedure is at least a factor 10 faster than the direct 
use of the true adjoint eigenfunction and is as accurate; it is 
a factor 2 more efficient than the difference computation of the 
group velocity and is much more accurate. 

(iv) Fast Matrix Setup -- A significant amount of computer, 
time in SALLY is spent setting up matrices to represent the 
Orr-Sommerf eld equation in terms of Chebyshev polynomial 
(about 40% of the computer time of the original SALLY code). 
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By careful reorganization of the code and new algorithms, we 
have been able to decrease this time substantially (to about 
20% of the run time of the new SALLY code) . 

Two new programs and the fast matrix multiply program 
discussed above are the essential ingredients of this speedup. 
A new Chebyshev derivative program was developed that, for 
typical LFC problems, runs about 20 times faster than the 
original code. It is based on the algorithm of Sec. 10 of Ref. 
Also, we have developed a new fast way to setup the matrices 
corresponding to non-uniform unperturbed flow. The new method 
is based on the convolutional structure of the equations. ^ 

(v) Spatial Eigenvalue Extrapolation — An improved 
technique to march from station to station across an LFC wing 
has been implemented. This new method is essential to get a 
good guess for the fast local eigenvalue solvers now used in 
SALLY. The idea is to use mean velocity profile information 
at two nearby stations to extrapolate the eigenvalue to the 
next station. This is done in a way that is transparent to 
the user of SALLY. 

(vi) Optimization Options — The original SALLY code 
offered the user only data on the most unstable temporal 
eigenvalue at a given frequency. The maximum is achieved over 
all wavelengths and propagation angles. Some new options have 
been added to the new code. They include: 

a) Fixed wavelength, fixed propagation angle. 

b) Fixed frequency, fixed wavelength. 

c) Fixed frequency, fixed propagation angle. 

In all cases, the user can determine the amplification of 
the given modes across the wing. 

(vii) Other Modi fications --In addition to the above modi- 
fications to the running SALLY code at NASA LaRC, a number 

of advanced developments have been made that are not yet part 
of the distributed version of SALLY. These include: spatial 



stability analysis; non-parallel flow analysis; and non-linear 
flow analysis. All these code modules are for three-dimensional 
boundary layers. 
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3 . Relation Between Temporal and Spatial Stability Theory 


We begin by reviewing the relation between temporal and 
spatial stability of a unidirectional, parallel, steady flow 
u(z). Let us study a three-dimensional mode that depends on 
x, y, t as exp(iax + i3y - iwt) . The Orr-Sommerfeld eigen- 
value problem leads to a relation between a, 3, U) of the form 

F(a,3,m) = 0 (3.1) 

where a, 3/ and a) are, in general, complex numbers and the 
function F is an analytic function of its arguments. In the 
following discussion, we use the notation that subscripts r 
and i indicate real and imaginary parts, respectively? for 
example, = Re (a). 

In the temporal theory, we set oc^ = 3^ = 0, so the mode 
is pure oscillatory in the streamwise and spanwise directions. 
Then (3.1) leads to two relations of the form 


= a) (a ,3 ) 

(3.2) 

r r r 

= ^ i <a r ,3 r ) 

(3.3) 


Transition prediction for LFC studies using an amplification 

factor criterion requires determination of the most unstable 

wave. To do this locally at each station on the wing, we maximize 

w . as a function of a and 3 . This maximum growth rate is 
l r r 

achieved when 


9 (jl) . 

i 

9a 


9oj . 


93 


- = 0 


(3.4) 


The two relations (3.4) determine the mode a r , 3 r of the most 

unstable three-dimensional temporal instability. 

The spatial stability theory for this unidirectional 

parallel flow is as follows. We set = 0 and solve (3.1) 

for a . and 3 • : 
i i 


6 



(3.5) 


“i = “i^r'W 

(a r , 3 r ^ r ) (3.6) 

The amplitude of this mode behaves like exp(-a^x - 3^y) , so 

the maximum growth rate is achieved by maximizing 
2 2 

a. +3. with respect to a , 3 , and to . However, the 
ii r r r 

Cauchy-Riemann equations imply that either = 0 or 3^ = 0. 

But a. — 0 is unrealistic. If a - = 0 then the wave is 

l l 

an 'edge 1 wave with its maximum growth perpendicular to the 
free stream. 

It follows that a reasonable procedure to treat spatial 
instability of a parallel flow is to choose = 3^ = 0 so 

that 


m = gj ( a , 3 ) 
r r r r 

(3.7) 

a i = “i<a r -e r ) 

(3.8) 


For maximum amplification, we must then determine a r and 

3 r by 

9a- 9a- 

n 


In this way, all modal parameters for the most unstable 
eigenfunction for a unidirectional parallel flow are determined. 

Now let us consider the changes in the above theory for 
general parallel flows in which the undisturbed flow is 
(u (z ,T) ,v (z ,T) , 0) in which we allow a slow time dependence 
on the long-time scale T. In this case, it is at first natural 
to use only temporal amplification. We assume a space-time 
dependence of the form 


exp(ia r x + i3 r y) exp (-i0 r (T) + (T) ) (3.10) 


where 
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(3.11) 


“r " WV T) ' 


03 . = 

1 


3 

"9T 


(T) 


The solution of the eigenvalue problem implies that 


0 iT = f(a r ' 3 r' T) (3.12) 

0 rT = 9<a r ,3 r ,T) (3.13) 

For any given wavenumbers a , 3 * Qj_ T = 0 and (3.12) define 

a neutral surface T = T, T (a ,3 )• The total amplification between 

N r r 

T._ and T is 
N 


cf(a r ,3 r /T) 


rT 

T„ ( ot 
N r 


3 ) 

r 


f (a r ,3 r ,T , )dT' 


(3.14) 


At any given T, the most unstable wave is determined by the 
equations 


3 o _ 3 g 

3o“ “ 3"3 
r r 


(3.15) 


Once a r and 3 r are found from (3.15) , (3.13) gives w . 

Next, let us consider a spatially varying unidirectional 
flow u(X,z), where X is a long-space scale in the flow 
direction. In this case, it is natural to study spatial 
amplification in which the space-time modal dependence is of 
the form 


exp(i0 r (X) - 0 i (X) ) exp (i3^y - iw t) (3.16) 

The Orr-Sommerfeld eigenvalue problem implies that 

-0. v = f (X, 3 , a)) (3.17) 

1 A IT 

0 r x = g(X,3 r ,u3) (3.18) 

In this case, a neutral curve can be defined by 0. = 0 and 
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(3.16) so that X = The total amplification from 

to X is given by 


o (X, 3 r »w r ) 


rX 

f (x' f 3 ,( 0 )dx' 

% 


(3.19) 


At any given X, the most unstable wave is determined by the 
equations 


9a _ 9a _ 
9 8 r 3 


(3.20) 


while the local wavenumber 0 v is then determined by (3.18) . 

rX 

Finally, let us consider a general spatially varying three- 
dimensional boundary layer in which the undisturbed velocity is 
(u(X,Y,z) ,v(X,Y,z) ,0) . In the original SALLY code, we considered 
temporally amplified modes on such a flow. Here we consider a 
self-consistent description of spatially amplified modes. In 
this case, we assume the space-time dependence of the mode to be 
of the form 


exp(i0 r (X,Y) - 0 i (X,Y)-iw r t) (3.21) 

The Orr-Sommerfe Id eigenvalue problem implies that 


iX f < X ' Y ' 0 rX ' G rY' aj r ) 

(3.22) 

iY = g ^ X ' Y ' 0 rX' 0 rY /OJ r ) 

(3.23) 


The problem now is that there is no completely rational 
way to decide which solution of (3.22) - (3.23) to choose for 
maximum amplification. In particular, there is no obvious 
reason why the direction of growth need coincide with either 
the free-stream flow direction or the normal to the modal 
wavefronts. It seems most reasonable to study instead a wave- 
packet in which the evolution of a given imposed perturbation 
near the leading edge of the wing is followed in its motion to 
the leading edge. This approach will be examined in more detail 
in later studies. Here we concentrate on the following simplified 
approach . 
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We assume that the direction of maximum growth -V0. is 
parallel to the freestream direction Ux + Vy. At a neutral 
surface, we must require that 0^ x = = 0; given the 

frequency oo r and wavenumbers © rX and 0 ^ , these equations 
determine a point X^, Y^. Using this fact, we can define 
various kinds of neutral curves. For example, those waves 
that initiate at angle $ with respect to the freesteam 
direction must satisfy 

Ua + VB 

—■ — ■ -- - - = cos <j> (3.24) 

/U 2 +V 2 /a 2 +3 2 
r r 


Eqs . (3.22) - (3.24) on the neutral surface allow elimination 

of a and 3 so we get a relation of the form 
r r 

H^XjY,^,^) = 0 (3.25) 


which is the equation of the neutral curve. 

The total amplification between Y N and X, Y is 

given by 


Q 


r (X,Y) 

Wv 


V0. -dX 

l 


(X , Y) 

(X N'V 



ds 


(3.26) 


which is independent of the path connecting X^, Y N with X, Y. 
Also, since -V©^ is assumed parallel to (U,V) , 


V 

U 


f (X,Y,a r ,3 r /W r ) 
g (X,Y,cT , 3 r / a> r ) 


(3.27) 


The maximum total amplification at a given point X, Y is 
achieved at a given frequency w as follows. Eq. (3.27) 
determines 3 r as a function of a r . Then starting from X, Y 
we march (3.26). and determine the most dangerous mode a by 

maximizing o . The integration in (3.26) stops as soon as the 
neutral curve (3.25) is reached. Once the maximum for a given 
frequency is found, the maximum amplification rate over all 
frequencies is found. 
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4. 


New Methods for Stability Analysis of Compressible Flows. 


The straightforward extension of the SALLY stability analysis 
code to include compressible effects can lead to highly inefficient 
use of computer resources. While the general three-dimensional 
disturbance of an incompressible parallel flow satisfies the 
fourth -order Orr- Somme rf eld equation, no single equation describing 
the evolution of even two-dimensional general compressible 
disturbances has yet been formulated. Since SALLY applies matrix 
methods (in order to be a 'black box' not requiring much user 
interaction) , the absence of a single stability equation for 

compressible flows implies that computer memory requirements 

2 3 

increase by K while computer time increases by K , where K 

is the number of equations entering the system that must be 

solved. Two-dimensional disturbances satisfy a sixth-order 

system that is easy to reduce to K = 3 second-order differential 

equations; three-dimensional disturbances satisfy an eighth-order 

system that is easy to reduce to K = 4 second -order equations. 

In both cases, the degration of computer time within the current 

SALLY code is unacceptable. 

We have developed a new method that allows easy solution 

of a very general class of stability problems by the spectral 

methods employed in SALLY while maintaining the accuracy and 

efficiency that spectral methods offer relative to more 

conventional difference methods. The key idea is a general 

algorithm to reduce an arbitrary system of K first-order 

'th 

differential equations into a single, explicit K~ — order 
differential equation. Since, in contrast to difference methods, 
the accuracy of spectral methods does not degrade seriously with 
increasing order, while the matrix methods employed by SALLY 
do degrade computationally with increasing number of equations, 
the reduction to a single high-order equation provides an 
attractive general way to solve the compressible flow stability 
problem efficiently. As a by-product of our analysis, we obtain 
a single equation for the evolution of general three-dimensional 
disturbances in a three-dimensional compressible boundary layer. 
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It is easy to construct a system of K first order 

til 

differential equations from the single K — order linear 
equation 


,K 


dx 


y 0 

IT + a K-l (x) 


,K-1 

d y, 


dx 


K 


+ + a Q (x)y 0 (x) = 0 


If we set y^ (x) 


d^y/dx k (k = 0,1 


,K-1) then 


(4,1) 


^=y 

dx Y k+1 


k = 0, . . . ,K-2 


(4.2) 


dy K-l 

___ = -a K _ 1 y K _ 1 - a K _ 2 Y K-2 - ' * *" a 0 Y 0 

Conversely, if the K functions y^ (k = 0,...,K-1) 
satisfy the general linear system 

(x) y (4.3) 


then it is usually true (and is true for the flow stability 

th 

analyses considered here) that yg (x) satisfies a K — order 
linear equation of the form (4.1). 4 To see this, we need 

only introduce the differentiated systems 

,P 

— p y = & (x)y (4.4) 

dx 1 


where 


dA (x) 

^P+l (x) = §p (x) § 1 (x) + — (P - 1) (4.5) 


Next, we solve (4.4) for y ^ ' Y 2 ' * * * ' Y K-l 3S -^ unct: *- ons °f 

Y 0 /Y 0 ,Y 0 ' ’ * ’ 'd K ^ and substitute into (4.4) with p = K 

to obtain a single K— order equation for y^ . 

There are some additional features that must be considered 
when the above algorithm is applied to compressible flow stability 
analysis. First, a mapping of the boundary layer direction must 
be made so the system of differential equations is a singular 
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system when expressed in the form (4.3). It is better to proceed 
as follows. Let A be the matrix of coefficients of the 
linearized compressible Navier-Stokes equations so that the 
boundary layer stability problem is of the form (4.3) with 
0 £ x < 00 together with certain boundary conditions at x = 0 
and x -»■ 00 . Then we make the algebraic mapping 

z = 2x/(x+L) - 1 (4.6) 

where L is a scale parameter and -1 <_ z < 1. Then (4.3) 
becomes 


ZDy = Ay (4.7) 

2 

where D = d/dz and Z = (1-z) /2L . Then the algorithm given 
above constructs a single K— order equation for y^ in terms 
of the operator ZD. We have developed an efficient program 
that converts this latter equation to an equation for y^ in 
terms of the differentiation operator D alone. Other new 
features include the automatic treatment of boundary conditions 
and the eigenvalue solver, to be discussed below. 

We have developed a computer program that solves the general 
problem (4.7) with rather arbitrary boundary conditions. Here 
we outline the organization of this code. A user-supplied sub- 
routine provides the coefficients of the matrix A in (4.7) . 

~ th. 

Then a subroutine calculates the coefficients f. of the K— 

k 

order equation 

(ZD) K y 0 + f K _ 1 (ZD) K_1 y 0 +...+f o y o = 0 (4.8) 

This subroutine also expresses y^,...,y K _^ in terms of y Q 
and its derivatives at the boundaries z = +1 so that arbitrary 
boundary conditions in terms of y n » • •• »y v can be applied. 

This subroutine is independent of the specific equations being 
solved. The output of this subroutine is a single matrix 
representing (in Chebyshev coordinates) the effect of the 
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differential operator (4.8) on y^ including both boundary 
conditions and the reduction of the operator ZD to D. 

The final step of the program is to compute the eigenvalues 
using the fast local eigenvalue subroutine already resident in 
SALLY. Since we use a local eigenvalue routine, it is necessary 
to provide it with an approximation to the eigenvalue. In a 
compressible flow problem this guess is best obtained either 
from a previous run or by slowly increasing the Mach number from 
an incompressible flow case in which a global eigenvalue routine 
like that employed by the incompressible SALLY code works 
efficiently . 

In summary, our new program provides a very general solution 
to eigenvalue problems of almost arbitrary complexity. For the 
three-dimensional modes of three-dimensional boundary layers, 
our new general program works only about 30% slower than the 
current SALLY code which can handle only incompressible flows. 
Another nice feature of the new program is that it is written 
in such a way that new problems can be solved almost trivially 
merely by inserting new subroutines for the coefficient matrix A. 
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5 . Spectral Methods for Boundary Layer Equations . 


In this Section, we present a brief description of some 

work we have done applying spectral methods to nonsimilar laminar 

5 . . 6 

boundary layers. The basic equation is 


3 3 f + f & + g (5) 


9n 


dr\2 


,9f.2l _ ^f9f 9 2 f 9 f 9 2 f 

W J 9?9n “ 8? a n 2_ 


(5.1) 


where f is the dimensionless streamf unction, £ and r\ are 
the Levy-Lees transformed variables, and the boundary conditions 
are 


f(?,0) = ff(S,0) = 0 §fU/D "■ 1 (riw) (5 - 2) 

We apply spectral-Chebyshev methods in n to the test case of 
Howarth 1 s flow in which 


S(5) = (5.3) 

To apply spectral methods in £, we first truncate 
0 < n < “ to 0 n R and then apply the linear mapping 

s = 2r\/R - 1 (5.4) 

to transform the range of the independent variable to -1 £ s £ 1. 
The boundary conditions become 

f (?> = !f(s-o) = °- H (5 ' R> = 1 (5 - 5) 

We use a spectral collocation method is s based on 
expansion of f in a series of Chebyshev polynomials T^(s); 
in the ^-direction, we use a Crank-Nicolson scheme like that 
employed by Keller and Cebeci. ^ 

In Tables 1-2, we compare our results with those of Keller 

f . : 

and Cebeci. It is apparent that we achieve limiting accuracy 
in the spectral direction r) with R = 6 and only 19 Chebyshev 
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TABLE 1. SPECTRAL METHODS FOR HOWARTH’S FLOW 
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for Howarth's Flow 


TABLE 2. Finite Difference Results 0 


Truncation R 

6 

# pts in r| 

19 

# pts in £ 

16 

Total // pts 

304 

A 

0.21063 

B 

0.41647 

Best Richardson 
extrapolation on 
this grid: # of 
pts used in 
extrapolant 

346 

A 

0.20735 

B 

0.03511 


6 

6 

61 

121 

16 

51 

976 

6171 

0.207232 

0.207227 

0.030932 

0.030832 


1708 

9282 

0.207410 

0.2071427 

0.031213 

0.0305307 
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polynomials. It is also apparent that the spectral method 
achieves results of at least the accuracy achieved by the 
extrapolation methods with an order-of-magnitude fewer grid 
points. In particular, the results obtained by Richardson 
extrapolation in the £ -direction of the spectral scheme uses 
only 253 points and is nearly as accurate as that obtained by 
Keller and Cebeci with 9282 points. The number of points used 
in the spectral calculation could be reduced by choosing R = 6 
and using only 19 polynomials (a total of 209 degrees of freedom) 
without any expected loss of accuracy. 

Spectral methods for boundary layer problems offer much 

promise. They seem to require an order-of-magnitude less 

computer memory than comparable finite-difference schemes and 

at least several times less computer time to achieve accurate 

results. Efficient methods to solve the spectral equations are 

5 

discussed by Orszag and Jayne. 
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6 . Nonlinear Nonparallel Stability Analysis 


In Sec. 2 of Ref. 2, the equations of two-dimensional 
nonlinear, nonparallel disturbances in an incompressible 
boundary layer are formulated. We have extended these equations 
to three-dimensions and have developed computer programs, con- 
sistent with SALLY to solve the resulting equations. 

Several simplifications in the equations of Ref. 2 have 
been made. As discussed in Sec. 2 of this report, it is not 
necessary to use the true adjoint eigenfunction of the Orr - 
Sommerfeld equation to impose solvability conditions, like that 
used to derive the group velocity. Instead, the matrix adjoint 
available from the subroutine EIGREV works just as well as 
is available with no extra computing effort. The run time for 
a typical three-dimensional nonlinear, nonparallel stability 
calculation is about twice the computer time of the linear 
version of the SALLY code. 

To date, we have run two sets of calculations with the 
full nonlinear, nonparallel version of SALLY. First, we have 
tested the formulation of the nonparallel flow terms by 
comparison with the results of Saric and Nayfeh.^ Agreement 
with the published data was achieved. Second, we have com- 
pared the results of SALLY with the results obtained by Antar 
and Collins^ for the nonlinear three-dimensional evolution 
of a Blasius boundary layer. Good numerical agreement was ob- 
tained. We have compared these results on the flat plate 
boundary layer with direct numerical simulations of transition 
in this flow. ^ The nonlinear stability theory calculations pre- 
dict a somewhat slower growth of disturbances than given by the 
full Navier-Stokes calculations. 

It is recommended that future runs with the SALLY code 
produce the nonlinear, nonparallel information available from 
the present extension of the code. In particular, it has been 
shown by Orszag^ and Orszag and Kells^ 0 that three-dimensional 
nonlinear effects are crucial in the transition process. 
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